arXiv:1501.06926vl [physics.soc-ph] 27 Jan 2015 


Published, in New J. Phys. 17, 015012 (2015) 


Comparative analysis of existing models for power-grid synchronization 

Takashi Nishikawa 1, a ) and Adilson E. Motter 1 

Department of Physics & Astronomy and Northwestern Institute on Complex Systems, 
Northwestern University, Evanston, IL 60208, USA 

The dynamics of power-grid networks is becoming an increasingly active area of re¬ 
search within the physics and network science communities. The results from such 
studies are typically insightful and illustrative, but are often based on simplifying 
assumptions that can be either difficult to assess or not fully justified for realistic ap¬ 
plications. Here we perform a comprehensive comparative analysis of three leading 
models recently used to study synchronization dynamics in power-grid networks—a 
fundamental problem of practical significance given that frequency synchronization 
of all power generators in the same interconnection is a necessary condition for a 
power grid to operate. We show that each of these models can be derived from first 
principles within a common framework based on the classical model of a generator, 
thereby clarifying all assumptions involved. This framework allows us to view power 
grids as complex networks of coupled second-order phase oscillators with both forc¬ 
ing and damping terms. Using simple illustrative examples, test systems, and real 
power-grid datasets, we study the inherent frequencies of the oscillators as well as 
their coupling structure, comparing across the different models. We demonstrate, in 
particular, that if the network structure is not homogeneous, generators with identical 
parameters need to be modeled as non-identical oscillators in general. We also discuss 
an approach to estimate the required (dynamical) parameters that are unavailable 
in typical power-grid datasets, their use for computing the constants of each of the 
three models, and an open-source MATLAB toolbox that we provide for these com¬ 
putations (available at https://sourceforge.net/projects/pg-sync-models). 
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I. INTRODUCTION 


In the field of network dynamics, power grids have long served as one of the main model sys¬ 
tems motivating the study of synchronization phenomena 1-3 . As the field became more mature, 
an increasing number of researchers have been seeking to apply ideas from past theoretical studies 
to power grid-specific problems 4-10 . The need for such applications comes from the fact that, de¬ 
spite the extensive engineering literature on power systems, there remains a largely under-explored 
problem of how the large-scale network structure influences the collective dynamics in power grids. 
While previous studies have emphasized the detailed modeling of relatively small test systems, the 
increasing availability of data processing tools, substantial computing power, and theoretical de¬ 
velopments in network synchronization are now making it possible to address large-scale properties 
of power-grid systems. 

A major concern for power grids is the stability of desired states, in particular synchronization 
stability of the power generators, which is a condition required for their normal operation. A 
frequency-synchronous state of n g power generators is characterized by 

Si = 62 = ■ ■ ■ = 6 ng , ( 1 ) 

where Si = Si(t) is the angle of rotation associated with the z-th generator at time t. Studying 
the stability of synchronous states of an alternating current interconnection against perturbations 
requires a network model capable of describing the coupled dynamics of power generators. Different 
models have been used in different publications in the physics literature, and there has been no 
comprehensive comparison to clarify how these models are related to each other. Providing such a 
comparative analysis is the main focus of this article. 

Here, we discuss three leading models, which we refer to as the effective network (EN) model 8,11 , 
the structure-preserving (SP) model 12 , and the synchronous motor (SM) model 13 . Each model is 
described as a network of coupled phase oscillators whose dynamics is governed by equations of 
the form 

9/t. r). 

—-Si H- -Si = At- V Kij sin(<Jj - Sj - 7 *_,•), ( 2 ) 

ur ur , 7 - U . 

where ujr is a reference angular frequency for the system, and Hi and Di are inertia and damp¬ 
ing constants characterizing the oscillators, respectively. The differences in the three models are 
reflected in the definition of the parameters Ai , , and 7 ^, as well as in the number of phase os¬ 
cillators. The phase angle Si of oscillator i may represent either a generator or load. While all three 
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FIG. 1. Modeling of power-grid network dynamics. (A) Simple network representation with nodes 
representing generators or loads, and links representing transmission lines or transformers. Here 
we used the example case consisting of two generators (nodes 1 and 2) and one load node (node 3), 
which is discussed in Section VI A. (B) Representation of the electrical properties of the components 
in the same network. There are three possible ways to represent the load node, which are used in 
three different models. (C) Representation of the system as a network of coupled oscillators for 
each choice of load representation in (B). For the system parameter values given in Section VIA, 
the network dynamics obeys Eq. (2) [or equivalently, Eqs. (28), (29), and (30) for the EN, SP, and 
SM model, respectively] with the indicated values of Aj, Kij, and 'fij. Each of the three dynamical 
models has its own definition of nodes that are different from that used in (A). In (B), these nodes 
are shown as black dots and indicated by orange, blue, and green indices for the EN, SP, and SM 
models, respectively. The same coloring scheme is used for the node indices in (C). Note that in the 
SP model the generator terminals are treated as load nodes with zero power consumption (nodes 
3 and 4), separately from the generator internal nodes (nodes 1 and 2), leading to the 5-node 
representation. 
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models represent the n g generators as oscillators, they are distinguished mainly by their different 
modeling approaches for the loads, which are representations of individual or aggregated consumers 
that draw power from specific points in the (transmission) network. The EN model represent the 
loads as constant impedances rather than oscillators, thus focusing on the synchronization of gen¬ 
erators as second-order oscillators. The SP model represents all load nodes as first-order oscillators 
(i.e., Hi = 0 ), and each generator is represented by two oscillators, including one for its terminal 
(a point connecting the generator to the rest of the network). The SM model assumes that the 
loads are synchronous motors that are represented as second-order oscillators. Figure 1 shows a 
simple example network that illustrates how these differences lead to different network dynamics 
representations with different number and values of model parameters. Note that the parameters 
are denoted in the figure using appropriate superscripts (EN, SP, and SM), which is a convention 
we will use throughout the article. 

The parameter A%, along with Di, determines the i-th oscillator’s inherent frequency, u* := 
wr( 1 + Ai/Df), which is the equilibrium frequency of oscillator i in the absence of the coupling 
term [the summation term in Eq. (2)]. For generators, this frequency is typically much higher than 
the system reference frequency ujr, as the two terms on the r.h.s. of Eq. (2) balance each other in a 
steady state. In a realistic setting, however, the instantaneous frequency of a generator is unlikely 
to actually reach this inherent frequency, since the system operator would take control actions 
well before the frequency deviates too far from the designated system frequency ojr. In addition, 
the equation of motion for generators and motors we derive below assumes that the frequency 
remains close to ujr and thus would no longer be valid if the frequency deviates too far from ujr. 
Nevertheless, this definition of inherent frequency can be useful in characterizing the nature of 
individual oscillators and, in fact, is analogous to the one usually used in studying networks of 
coupled phase oscillators. Note that “natural frequency” is a term commonly used to refer to u ;* in 
that context, but it has a very different meaning in the power systems literature; see Appendix A. 
The parameter K t] > 0 represents the strength of dynamical coupling between oscillators i and j, 
while 7 ij represents a phase shift involved in this coupling. We note that Eq. (2) can be regarded as 
a second-order analog of the Kuramoto model 14 with arbitrary coupling structure. Similar forms of 
second-order equations for coupled oscillators have been used to study synchronization and phase 
transitions outside the context of power grids 15-18 . 

The paper is organized as follows. In Sec. II, we discuss how to characterize and compute 
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the steady state of a power-grid network. We then turn our attention in Sec. Ill to the problem 
of describing the dynamics of individual generators, deriving the basic equation of motion and 
describing the electric circuit representation of a generator under the classical model. In Sec IV, we 
derive the EN, SP, and SM models and discuss their differences and similarities. In Section V, we 
describe how to obtain model parameters required for numerical simulations. Finally, we discuss in 
Sec. VI an instructive example system and a selection of test systems and real power-grid datasets, 
focusing particularly on the heterogeneity of the constants Ai and the inherent frequencies as well 
as on the coupling structure encoded in the constants Kij. We make concluding remarks in Sec. VII. 


II. STEADY STATES OF POWER-GRID NETWORKS: THE POWER FLOW 
EQUATIONS 

A power grid is a system of electrically coupled devices whose purpose is to deliver power from 
generators to consumers. To view the system as a network, we define a node to be a point in the 
system at which power is injected by a generator or extracted by power consumers, or a branching 
point through which power is redistributed among the transmission lines connected to the point. 
A link is then defined as an electrical connection between a pair of such nodes and can represent 
a transmission line or a transformer. This network representation of the physical structure of 
a power grid is illustrated in Fig. 1A using a simple example of two power-injecting nodes (1 
and 2) connected to a single power-consuming node (3) by two transmission lines. Transmission 
lines have impedances, which are represented by complex numbers in general, and their parallel 
conductors have capacitance between them. In power system analysis, a transmission line is usually 
represented by the so-called n model, in which the two nodes are connected by an impedance, 
with two capacitors (of equal capacitance) connecting both sides of the impedance to the ground. 
This model is illustrated for the transmission line 1-3 in Fig. IB, which has two capacitors with 
impedance 1 /(j&13/2), 613 > 0, where we denote the imaginary unit by j := \J— 1. A transformer is 
represented by a model in which the (complex-valued) voltages on the two sides of the transformer 
maintain a constant ratio (generally complex-valued to account for a possible voltage phase shift). 
The standard approach 19 , which we adopt here, is to represent the parameters of these models for 
transmission lines and transformers in terms of equivalent admittances (the inverse of impedances). 
The structure of the entire physical network can then be represented 19 by the (complex-valued) 
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admittance matrix Yo = (Yo^), where an off-diagonal element Yoij is the negative of the admittance 
between nodes i and j i and a diagonal element Yqh is the sum of all admittances connected 
to node i (including the shunt admittances to the ground, which are parts of the models for 
transmission lines and transformers). 

The operating condition of a power grid can be characterized by the distribution of power flow 
through the physical network of transmission lines and transformers. In an alternating current 
grid, power is represented as a complex number, whose real and imaginary components, P, and 
Qi, are called active and reactive power, respectively. Alternatively, the power flow state can be 
characterized uniquely by the complex voltage V) (with magnitude | Vj| and phase angle fa, so that 
Vi = iVije- 7 ^*) and the complex power P r + jQi injected into each node i in the network. Given the 
locations and parameters of power generators and loads, as well as the parameters of the network 
summarized in the admittance matrix Yo, fundamental laws of physics—the Kirchhoff’s laws— 
determine the power flow state of the system. The laws can be used to derive the so-called power 
flow equations 19 , which provide the foundation for steady-state power system analysis: 

n 

Pi = 'Yh lWo«l sin (0j - - 7o ij), i = l,...,n, (3) 

j = 1 

n 

Qi = \ViVjY 0ij \cos(fa - fa - yo ij), i = l,...,n, (4) 

3 = 1 

where the complex admittances are represented in polar form as Yoij = I Yoij I e 3a ° i:i , we define 
7 o ij := oiQij — vr/2, and n is the number of nodes. To solve this set of 2 n nonlinear equations for all 
4n quantities that determine the power flow state of the system, we require that two real quantities 
per node are given as parameters. The usual assumptions are: 

• If node i is a generator node, the values of Pi and | V | are given as Pi = P* } i and U: = | Vj* |. 
This is reasonable because in real power grids the generators are scheduled to produce 
constant active power at a given time, and constant voltage magnitude is usually maintained 
by voltage regulators. Usually, one generator node is chosen to be a reference node, for which 
i fa is set to zero instead of specifying the value of P,;. For the steady-state analysis, we need 
at least one node with unspecified P; because the total power generation cannot be known a 
priori even when the total power consumption is known, since the difference, which equals 
the power lost in the transmission lines and transformers, depends on the power flow state. 

• If node i is a load node, the values of Pi and Qi are given as P* = —P/j and Qi = —Q\i- 
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For positive active and reactive power consumptions, Pi and Qi take negative values. These 
values can be determined from historical data, load forecast, or measurement 11,19 . 

Given the parameter values and the admittance matrix Yo, the nonlinear equations (3) and (4) are 
solved numerically, which can be performed by a variety of software, including MATLAB-based, 
freely available packages such as Power System Toolbox (PST) 20 and MATPOWER 21 . In the 
following, we will denote the values of the state variables for a power flow solution as P*, Q*,V*, 
and </>*. 

III. DYNAMICS OF INDIVIDUAL GENERATORS WITHIN THE NETWORK 

A. Mechanical representation — The swing equation 

Many components of a power grid, such as generators and loads, are dynamic and the state 
of the component can change in time. A dynamical model of the grid is obtained when the two 
parameters given for each node to determine the power flow state are replaced by a set of differential 
equations that, together with the power flow equations [Eqs. (3) and (4)], allow the determination 
of the state of the system as a function of time. Instead of the constant power injection and voltage 
magnitude at a generator node, one can write an equation of motion describing the dynamics of 
the generator rotor. 

Most generators in today’s power grids have a rotating mass (a rotor) that is driven by mechan¬ 
ical torque to generate electrical power. There exist models with various levels of sophistication 
and complexity, but in all such models, the rotational dynamics of the rotor is ultimately governed 
by a fundamental law of physics: the Newton’s second law. The electrical output from such a 
generator is coupled to other generators in the grid through a network of transmission lines, trans¬ 
formers, and other devices, which serve to transport the generated power to consumers. The power 
demanded by the network is felt by the generator’s rotor as an electrical torque, which is usually 
a decelerating torque. It is the balance between the mechanical input power and the electrical 
output power that determines the dynamics of the rotor. 

To derive the equation of motion governing the dynamics of such a generator, we set the rate 
of change of the angular momentum of the rotor equal to the net torque acting on the rotor: 

JS = T m - D m u> - - D e Aw - T e , (5) 
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where J is the moment of inertia in kg-m 2 , 5 is the angle of the rotor relative to a frame rotating at 
the reference frequency ujr in rad/s, T rn is the mechanical torque in N-m driving the rotor, D m is 
the damping coefficient in N-m-s for mechanical friction, u is the angular frequency of the rotor in 
rad/s, R is the regulation parameter in rad/N-m-s characterizing the proportional frequency control 
by a governor, A uj := oj — ujr is the frequency deviation, D e is the damping coefficient in N-m-s for 
the electrical effect of the generator’s damper windings, and T e is the typically decelerating torque 
in N-m due to electrical load in the network. Noting that 5 = u> — ur = Aw, we can rewrite the 
equation as 

JS + D5 = T m - T e , (6) 

where D := D m + D e + 1/R and T m := T m — D m ojR is the net mechanical torque, accounting for 
frictional loss at the reference frequency. Multiplying both sides by u and using the fact that torque 
in N-m multiplied by angular velocity in radians per second gives power in watts, the equation can 
be written in terms of power: 

Jlur5 + Duj r 5 = ^ (T m u - T e uj) « P m - P e , (7) 

UJ 

where we define P m := T m uj and P e := T e u, and we assumed the factor ojr/oj to be nearly equal 
to one, i.e., that the generator’s frequency u is close to the reference frequency ur. 

This approximation, which can be formalized using singular perturbation analysis 22 , is valid 
and considered appropriate for studying power system stability, where the key question is whether 
the system desynchronizes after a perturbation from a steady-state operation near the reference 
frequency. If, for example, a disturbance leads to even 1% deviation of a generator’s frequency 
from the 60 Hz reference frequency for a period of just one second, it would lead to an increase 
in angle difference equivalent to more than half a rotation. This much of angle deviation is highly 
likely to cause loss of synchronization in practical situations 11 . While the approximation co ~ cor 
is widely accepted in this context, a more detailed electromechanical model that does not require 
this approximation is also available 2 ' 1 . 

We now divide both sides of the equation by the rated power Pr (used as a reference) to make 
P m and P e per unit (p.u.) quantities, which is a normalization procedure commonly used in power 
systems studies. The factor Jur then becomes 2H/ujr, where we defined the inertia constant H := 
W/Pr (in seconds) and the kinetic energy of the rotor W := \Ju 2 R (in joules). The factor Dujr 
becomes D/ur, where we defined the combined damping coefficient as D := Dujr/Pr (in seconds). 


We then obtain what is known in the power systems literature as the swing equation 11,19,22 : 

‘ITT D 

—6 + —5 = P m - P e , (8) 

LOR LOR 

which we use here as the fundamental equation of motion for a generator. The term P m represents 
the net mechanical power input to the rotor, while P e represents the electrical power demanded 
by the rest of the network and includes terms that depend explicitly on 6 and the state variables 
of the other generators and loads in the network. 


B. Electric circuit representation — The classical model 

In general, both P m and P e have nonlinear dynamics, which may need to be accounted for 
in high-accuracy simulations that power system engineers require. There has been substantial 
effort in the power systems literature to model the dynamical behavior of the generator’s internal 
magnetic flux, which affects the electrical power output Pe 11,19,23-26 . One may also need to include 
the nonlinear dynamics of the governor that controls the generator frequency and the excitation 
system that controls the voltage magnitude at the generator terminal. Here we aim to study power 
grids as a complex dynamical system and focus on the way in which the network structure influences 
the synchronization dynamics of generators in simplest settings that are nevertheless realistic. For 
this purpose, we now describe a model of a generator used commonly in the engineering literature, 
particularly in theoretical studies, and forms a basis of all three network models we present below. 
In this so-called classical model, a generator is represented as a voltage source with constant voltage 
magnitude \E\ connected to the terminal node through a reactance x' d > 0 (see Fig. IB). The phase 
angle of the voltage source is assumed to be equal to the rotor’s rotational angle and thus denoted 
by 5. In addition, the mechanical power input P m to the rotor in Eq. (8) is assumed to be constant. 
Constant voltage magnitude and mechanical power are approximations that are valid for short-term 
dynamics. For analysis of transient stability after a disturbance, the approximation is considered 
valid for simulation of the “first swing” of the resulting oscillation of uo or 5, which typically covers 
a time period on the order of one second 11 . The classical model sometimes refers to the one in 
which damping is ignored [i.e., D = 0 in Eq. (8)], but here we consider damping effects explicitly, 
as they can have non-negligible effect on the stability of steady-state power grid operation and can 
potentially be used as tunable parameters for optimizing the stability 8 . 
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In a real power grid under a steady-state condition, the system frequency is continuously mon¬ 
itored and is in fact actively controlled (by adjusting P m of its generators) to stay close to the 
reference frequency of the system. Because of this, in power system stability analysis, it is usually 
assumed that all generators are initially synchronized at the reference frequency ujr in a steady 
state. Then, for a given steady state in which the generator is delivering active power P* through 
its terminal, we have 5 = 5 = 0 and P m = P e = P*. This value of P m is then held constant when 
studying stability using the swing equation and the classical model of the generator. 

Regardless of whether the grid is under steady-state condition or not, the expression for the 
electrical output power is given by the so-called power-angle equation, P e = ^ E ^ sin(<5 — </>), where 
the complex voltage V = (Rle- 7 ^ at the terminal has generally time-dependent magnitude and 
phase, and E = \E*\e^ is the voltage at the internal node defined to be a point between the 
constant-magnitude voltage source (with time-dependent phase) and the transient reactance (see 
Fig. IB). Equation (8) then becomes 


^5+—5 = P*-t E ^ S m(5-cl ) ), (9) 

UR UR y X' d 

While H, D, and x' d are constants that characterize physical and electrical properties of the specific 
generator, P* and \E*\ are parameters that depend on the steady-state distribution of power flow 
across the network. By eliminating the current I from the expression for the complex power, 
P* + jQ* = VI (where the bar denotes the complex conjugate), and the Ohm’s law, jx' d I = 
E — V, applied to the transient reactance, the steady-state internal voltage magnitude \E*\ can be 
calculated as 


| E*\ 2 


P >'d 

W*\ 


+ 



( 10 ) 


where Q* is the reactive power injected by the generator into the terminal and |R*| is the terminal 
voltage magnitude in the steady state. We thus see that the state of the terminal node given by 
P*. Q* and |R*| determines some of the parameters required to simulate generator dynamics using 
Eq. (9). This is a crucial point that, as we will show below, has significant implications for the 
modeling of a power grid as a network of coupled oscillators. In a typical stability analysis using 
the classical model, P* and |R*| are given as input data, and Q* is obtained by solving the power 
flow equations (3) and (4) for the entire network, given necessary input data at other nodes. The 
power flow solution also provides values of the phase angles, 5* and </>*, in the steady state, which 
can be used as the initial condition for Eq. (9). Note, however, that the equation is not closed, as 
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the dynamics of 4> = 4>(t) and V = V (t) are not yet specified. It is mainly the phase angle <j> that 
serves as the medium through which this generator is coupled to the rest of the system. 


IV. COLLECTIVE DYNAMICS OF GENERATORS IN POWER-GRID 
NETWORKS 

The dynamical state of the network is characterized by four functions of time for each node: Pi, 
Qi, | Vi j, and </>$. The power flow equations (3) and (4), being valid at each instant of time for time- 
varying variables, provide two out of the four required equations per node for uniquely determining 
the dynamical state of the system. We thus need two additional equations/assumptions at each 
node. For a generator node, the swing equation effectively provides one equation, and the usual 
practice is to make the additional assumption that the reactive power Q g injected by the generator 
at its terminal node is constant over short time scales and equals its steady-state value Q*, which is 
computed from the (static) power flow equations. We would then have a complete set of equations 
that allows us to determine all four state variables of the generator, \E\, 6 , P g , and Q g , as a function 
of time, given the initial condition and the state variables (for all t) at all the other nodes. Note 
that E = \E\e^ is related to the terminal voltage V = \V\e^ through the transient reactance at 
all times. Of the four generator variables, 5 directly relates to Eq. (1) and thus is the most relevant 
for the problem of synchronization stability. 

Supplying the full set of equations for load nodes requires modeling of static and/or dynamic 
behavior of the loads. Load modeling is a hard problem because the power consumption at a node 
in a transmission network is typically an aggregate of a large number of loads of a wide variety 
of types and sizes, each of which can be time dependent and influenced by human activity. As 
a consequence, a number of different models have been proposed and used in the power system 
literature. As mentioned above, it is desirable to have a simple but realistic model that can clarify 
the role of network structure in the problem of power-grid synchronization stability. We thus 
provide in this article a comparative analysis of three approaches that can be used to recast the 
problem as that of synchronization in complex networks of coupled oscillators. These approaches 
lead to the three different models, the EN, SP, and SM models, which are all described by Eq. (2) 
but with different values and interpretations for the parameters Aj, Kij, and 7 ^. 
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A. The effective network model 


The dynamical interaction between a pair of generators is mediated by the paths of transmission 
lines connecting the generators and is expressed as the sinusoidal dependence on the terminal 
voltage phase in Eq. (9). While the nature of the coupling between generators is ultimately shaped 
by the structure of the transmission network and the distribution of loads, it would be insightful and 
useful if one could represent the coupling in a single term that depends on the state variables of the 
generators. This can indeed be achieved by modeling loads as constant impedances and reducing 
the physical network to what we call the effective network of interactions between generators. Since 
the model is derived through this reduction process, it is also called the network-reduction model 
or network-reduced model in the literature. 

In order to see how the reduction process leads to the effective network, which then determines 
the coupling constants Kij governing the network dynamics in Eq. (2), let n g and ni denote the 
number of generator terminal nodes and load nodes, respectively, so that n = n g +nn. In addition to 
these n nodes, we define an additional node to be a point between the internal transient reactance 
and the constant voltage source for each generator (see the two generators depicted in Fig. IB), 
making the total number of nodes N := 2 n g + ri£ = n g + n. By suitable re-indexing, we may assume 
that i = 1 , ... ,n g correspond to the generator internal nodes, i = n g + 1 ,..., 2 n g to generator 
terminal nodes, and i = 2 n g + 1 ,... ,N to load nodes. Define Y^ as the n g x n g diagonal matrix 
whose diagonal elements are the admittances (jx' d i) -1 ,..., {jx' dng )~ l of the generator transient 
reactances. We write the admittance matrix Yq for the physical network as 



where Yq 9 , Yq^, Yq 9 , and Yq 1 are the four blocks resulting from separating the first n g rows 
(columns) from the last rows (columns). 

The active power Pf t and reactive power Q\ i consumed at load node i (and possibly also at 
a generator terminal node) in a steady state is represented by an equivalent impedance to the 
ground 11 whose admittance is Yf :i = {Pf t — jQ\ j)/|V)*| 2 , where \V*\ is the voltage magnitude at 
the node. The constants Pf i and Q\ i are part of the input data for Eqs. (3) and (4) to solve for 
a power flow solution, which provides \V*\. The constant admittance Y),i is computed from these 
steady-state values and added to the corresponding diagonal components of Yq 9 and Yq to obtain 
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Y 99 and Yff. This representation requires the assumption that the power demand is constant. We 
thus consider the dynamics of the system on time scales that are short enough for the validity of 
the assumption, but long enough to address the problem of synchronization. We can now define an 
N x N admittance matrix Yq = (Y^T) that includes the links representing the transient reactances 
connecting the generators’ internal and terminal nodes, as well as the equivalent impedances for 
the loads: 



/ 

Y d 

—Y d 

0 

Y' ■- 

1 0 • — 


-Y d 

Y 99 +Y d 

Y 9 ^ 

1 0 


V 

0 

■y^S 

1 0 

<rU 

r 0 


( 12 ) 


where 0 denotes matrices (of appropriate sizes) whose elements are all zeros. 

Now let V 5 , V 1 , and be the vectors of node voltages for generator internal nodes, generator 
terminal nodes, and load nodes, respectively, and stack them vertically in that order to form the 
voltage vector V. Also let I 9 be the vector of currents injected into the system at the generator 
internal nodes. Because the loads are modeled as constant impedances to the ground, all nodes 
have zero injection currents except for the generator internal nodes. The Kirchhoff’s current law 
can then be written in the form I = YqV, or equivalently 


H 


Yd 

—Yd 

0 N 


H 


0 

= 

-Yd 

Y^ 9 + Y d 

1 0 



(13) 

\v 


l ° 

-y^ 9 

1 0 

Y u 

Y ° / 


w 



The system is then converted to I 9 = Y V 9 by eliminating V* and V*, a procedure known as Kron 
reduction 19 , where the resulting effective admittance matrix Y EN = (Y( EN ) is defined by 


Y en := Y'(l + Y^Y')" 1 , where Y' := Y 9 0 9 


Y 9i (Y 
1 o 1 1 1 


i — l yb? 

I 1 o ! 


(14) 


where 1 denotes the n g x n g identity matrix. The inverse of Y^ always exists, while the matrix 
(1 + Y ( ^ 1 Y , )~ 1 exists when x' di are small enough. The existence of the matrix (Yq 1 )^ 1 follows from 
the assumed uniqueness of the voltage vectors (with respect to a reference voltage). The symmetric 
n g x n g matrix Y EN defines an electrically equivalent network in which generator internal nodes 
i and j f i are connected by an effective admittance —Note that this notion of effective 
admittance is distinct from the inverse of the (two-point) effective impedance used in AC circuit 
theory 27 , but it is more suitable for our purpose here since it captures the idea that dynamical 
interactions between generators are determined by an effective network of admittances connecting 
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the generators. Also note that a similar effective network can be derived in the case of constant 
current injections at load nodes and generator terminal nodes [nonzero current vectors I 4 and 1 ^ 
replacing the zero entries on the left hand side of Eq. (13)]. 

Accounting for the transient reactances x' d ■ is important, since they are typically not small 19 
for real generators producing small power and vary widely from generator to generator. This is 
because the values need to be expressed in p.u. with respect to the system base (a common set of 
reference units for the system) when writing the admittance matrices Y^ and Yq. If instead the 
values were expressed in p.u. with respect to the rated power Pr of the generator, they tend to 
lie in a relatively narrow range 19 . A consequence of having x' d ■ >0 for all i is that the effective 
network represented by Y EN , which is obtained after eliminating all nodes except for the generator 
internal nodes, has the topology of a complete graph (Yjj 7 ^ 0 for all i and j). This follows from 
a general property of Kron reduction that two nodes are connected in the reduced network if and 
only if the two nodes are connected in the original network by a path in which all intermediate 
nodes are eliminated by the reduction process 28 . If one neglects the transient reactances and sets 
x'j • = 0, we see from Eq. (14) that we would have Y EN = Y'. The matrix Y 7 can be shown to 
equal the admittance matrix between the generator terminal nodes obtained by eliminating the 
load nodes through Kron reduction, which may or may not have the topology of a complete graph 
depending on the location of the load nodes. 

Using the effective admittance matrix Y EN , the single sinusoidal coupling term in Eq. (9) can 
be replaced by an expression for P e that comes from a power balance equation equivalent to Eq. (3) 
with \Vi\ replaced by \E*\. Yo by Y EN , and fa by Sf 

n g 

Pe,i = E \E*E*Y™\cos(5j -6i + af), (15) 

3 = 1 

where Yj EN = | Yj EN | exp(ja EN ) and all quantities must be expressed in p.u. with respect to the 
system base. This can be used to show that writing the swing equation ( 8 ) for each generator leads 
to an equation of the same form as Eq. (2), and the EN model is thus given by 
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(16) 


where G EN is the real part of the complex admittance Yf N - Notice that the number of phase 
oscillators in Eq. (2) in this case is n g , since generators are the only dynamical elements over the 
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relevant time scales and each generator is described by a single variable 8 % . The internal voltage 
magnitude \E*\ is a constant computed from a given steady-state power flow by applying Eq. (10) 
to each generator. Based on the same reasoning we used at the single-generator level in Sec. IllB, 
the constant net mechanical power input to the generators is equal to the steady-state electrical 
power output P* % . It can be seen from this equation that the parameter A EN [which determines the 
inherent frequency uj* = wr (1 + A EN /Dj)] and the coupling strength AT EN depend not only on the 
structure of the transmission network expressed in Yo, but also on the steady-state power flow over 
the network through the values of P*P^, Q} i , and V *, which affect Y EN and | E*\. The effective 
admittances in Y EN for realistic systems have imaginary parts that are mostly positive and much 
larger than the real parts 8 , indicating that inductive reactances are the dominant components of 
the effective network. This leads to phase shifts y EN ~ 0, and we see that the individual coupling 
terms in Eq. (2) tend to keep the angle differences between generators small, reflecting an inherent 
tendency of connected generators to synchronize with each other. This justifies the use of the 
assumption that all transmission lines are lossless (i.e., their impedances have no real components) 
and inductive (i.e., their impedances have positive imaginary components), since it implies that 
Re(l) EN ) = 0, Im(Y i EN ) > 0 and y EN = 0 for all i / j. Note, however, that this assumption does 
not imply that the diagonal components 1( EN have zero real components (G EN = 0), and thus does 
not imply that Af N = P* t . 

In a recent publication 8 , we derived a master stability function 29 for the EN model, which 
reduces the problem of synchronization stability against small perturbations to a simple spec¬ 
tral condition that involves both the effective network and the synchronous state. Based on this 
condition, we identified a systematic adjustment of generator parameters that enhances the syn¬ 
chronization stability and demonstrated the effectiveness of this approach using a selection of test 
systems and real power-grid datasets. Sufficient conditions for synchronization stability against 
larger perturbations have also been derived for the EN model, one in terms of the maximum 
voltage phase difference in an approximation of the corresponding power flow solution' (assuming 
lossless transmission lines), and another in terms of the minimum combined coupling strength (ac¬ 
counting for K t j , 7 ij, and Dj) among all node pairs 30 . A similar sufficient condition has also been 
derived for an extension of the EN model that incorporates a second-order model for the frequency 
control of individual generators 31 . It has been found 9 that certain structural motifs can cause the 
loss of synchronization stability in a special case of the EN model which assumes that all nodes 
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have generators with identical parameters and negligible x' di , all transmission lines are lossless and 
have equal reactance, and half of the generators (randomly chosen) have A EN = P and the other 
half have Af N = —P for some constant P > 0. Under these assumptions, it can be shown that 
Y en = Y' has the same topology as the network of transmission lines. This simplified version 
of the EN model corresponds to a power grid that is homogeneous in all aspects, except that the 
locations of power generators and consumers can be heterogeneous and the network topology can 
be arbitrary. 


B. The structure-preserving model 

A different approach is to seek to describe the dynamic behavior of the loads, which can lead 
to a model of the system that retains the structure of the physical network of transmission lines 
connecting the load nodes. In a steady state, the active and reactive power injected into load node i 
are constant: Pjy = P^ i and Qt d = Q\ i . The equivalent admittance we used for the loads in the EN 
model above means that the dynamic behavior of the power consumption by the load is described by 
a function of the time-dependent node voltage magnitude: P^ = G^i\Vi(t)\ 2 and Qi^ = -B^j|Uj(f)| 2 , 
where Ga and B^ are the real and imaginary parts of the constant admittance Yg t . In general, 
the power consumption can depend nonlinearly on both the node voltage magnitude |V)(i)| and 
phase frequency 4>i(t). The SP model uses an alternative model for the dynamic behavior of active 
power consumption, Pf d = PJ t + ^<5j(t), where Di > 0 is a constant and 5i(t) := is 

the voltage phase relative to a common reference frame rotating at the synchronous frequency 12 . 
This model results from the linearization of the power-frequency relation, assuming constant node 
voltage magnitude |V)| = |V)*| and small deviation 8 t from the steady-state frequency. The reactive 
power consumption is assumed to be constant, Q?^ = Q* ti . This load model leads to the equation 
of the same form as Eq. (2), where we set Hi = 0 and re-interpret Di as the linear coefficient in the 
frequency dependence of the active power consumption. The structure of the network connecting 
the generators and loads is given by an N x N matrix similar to Yq in Eq. (12): 
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The equation for each generator is also of the same form but with Hi > 0, and is in fact exactly 
Eq. (9) because the only link from the generator’s internal node is to its terminal node through a 
purely imaginary admittance, Y^^_ Ug = —(jx' d J _1 , i = 1,... , n 5 . Putting together, the SP model 
is given by 

2 M D • 

—-5 + —6 = Af p - KfJ + sin- S i+ng ), i = 1,... ,n g , 

UR Ur 

4SP p* j-SP \F*V*/A I 

• g,i'> XY i,i+n g * I v i / 


for generator internal nodes and 

N 

Y Kff sin(<5 ? ; - Sj - 7® P ), i = n g + 1,..., IV, 


—= 4 P - 

UR 


j=n g +l,j^i 
SP ._ zd* it/-*i2/^ySP 7^SP 


— _ P* _ — IT/*T/*V.^ P I i — ? _ n i r — n — n 


(19) 


7*7 :=a if - ^ Y if = l*ijH ex P(j«f/) 


for load nodes (including the generator terminal nodes), where Gf p denotes the real part of . 
Note that these equations have the same form as Eq. (2), and in this case the number of phase 
oscillators is IV = n g + n, since each of the generator internal and terminal nodes, as well as load 
nodes, is represented by a single variable <7 In contrast to the EN model, the SP model has the 
advantage that the physical structure of the transmission network represented by Yo is preserved 
as part of the admittance matrix Y SP and reflected in the coupling constants Kf p . This, however, 
comes at the expense of additional complexity associated with the larger set of equations (real 
power grids can have many more nodes than the number of generators) and increased uncertainty 
associated with the additional parameters Dj that must be estimated. Note that in the SP model 
the transient reactances x d - are part of the preserved network structure. Since x' d i are not negligible 
in general (as we discussed in Sec. IV A for the EN model), the corresponding links in the preserved 
network structure cannot be ignored, even though they do not represent any physical connections. 
The inherent frequencies uj* = ur( 1 + Af p /D{) are well defined in this model for both the internal 
and terminal nodes of the generators as well as for the load nodes. 

If all transmission lines are assumed to be lossless and inductive (which is well-justified), the 
admittance matrix Y SP is purely imaginary with positive imaginary component, which implies that 
7 ^ p = 0 for all i 7 j- Under this assumption, the SP model has been used 4 to design incremental 
adjustments to the active power output of generators or transmission line reactances to improve the 
linear stability of synchronous states. The sufficient synchronization condition established in Ref. 7 
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and mentioned in Sec. IV A for the EN model (also assuming lossless transmission lines) has been 
shown in the same publication to be applicable also to the SP model. Variations of the SP model 
that represent droop inverters as first-order oscillators 32 and incorporate stochastic deviations of 
frequencies 33 have been used to study synchronization stability in networks with renewable energy 
sources. 


C. The synchronous motor model 


Another way to model the dynamics of the loads while preserving the physical network structure 
is to use synchronous motors to represent the loads. A synchronous motor is essentially the same 
type of machine as a generator, except that the flow of power is in the opposite direction, converting 
electrical power into mechanical power, and thus can be modeled by the same swing equation, 
Eq. (8) with P m < 0 and P e < 0. In addition to the n g internal nodes we defined for the generators 
in the EN model, in this case one defines ri£ internal nodes for the synchronous motors representing 
the loads, which makes the total number of nodes 2 n. Since the motors are modeled in exactly the 
same way as the generators, this network model is mathematically equivalent to the EN model in 
which all n nodes of the physical network has “generators,” ri£ of which have negative mechanical 
power P m ,2 ^ 0. The matrix Yq will be replaced by 


V" — 
In- — 


Y' 
1 d 


-Y 7 

1 d 


-Y' 


Yn + Y' 


( 20 ) 


in this case, where Y^ is the n x n diagonal matrix whose diagonal elements are the admittances 
(J x 'd i) -1 ; • • •, (j x d n )~ 1 of the transient reactances of both the generators and motors. Kron re¬ 
duction is necessary to eliminate the terminal nodes of all generators and motors, from which we 
obtain an n x n matrix Y SM and a system of n coupled phase oscillators in the same form as 
Eq. (2). Note that, for a given system, the matrix Y SM for this model is in general different from 
that for the EN model. The generator/motor transient reactances are typically non-negligible, and 
the network represented by Y SM can be shown to have the topology of a complete graph when 
x'd j > 0 for all generators and motors, applying the same argument we used for the EN model. 
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The equations of motion for the SM model read 


—Si + —Si = Af M - Af M sm($i - Sj - 7g M ), i = l,...,n, 

Wfl UR . 


( 21 ) 
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where here A® M := P* ?; — |A*| 2 G Y | M for generators (i = 1,n g ) and A? M := ~P( % — |A*| 2 G| m 
for loads/motors [i = n g + 1,..., n). The internal node voltage magnitude \E*\ for the generators 
and motors is determined in the same way as in the EN model. While the representation of the 
loads as synchronous motors allows for a simplified description of the system with second-order 
differential equations for both generators and loads (in contrast to the SP model, which has both 
first- and second-order equations), such representation is typically not adopted in the power system 
engineering literature. Most of motor loads in the US are indeed induction motors, not synchronous 
motors 34 . 

As in the case of the SP model, the well-justified assumption of lossless and inductive trans¬ 
mission lines implies that = 0 for all i ^ j. This assumption is adopted in most of the 
previous studies that use the SM model, including the one that originally proposed the model 13 . 
The SM model has been used to show that decentralization of power generation makes synchroniza¬ 
tion more sensitive to temporary increase in power demand, while simultaneously making it more 
robust against removal of single transmission lines 5 . While transmission line removals naturally 
destabilize power-grid synchronization, it has also been shown using the same model that adding 
transmission lines can sometimes have the same effect 35 . The impact of decentralization has also 
been shown 36 to increase the order parameter that measures the degree of synchronization for a 
range of different random network topologies. In conjunction with the network topology of the 
European power grids, the SM model has been used to study the minimum coupling strength re¬ 
quired for synchronization 6 . Reference 37 studies synchronization in Kuramoto oscillator networks 
with bipolar distribution of inherent frequencies as a model of power-grid networks (equivalent to 
assuming Hi = 0 in the SM model), focusing on the effect of correlation between the inherent 
frequencies of two nodes that are connected. All of the above studies make additional simplifying 
assumptions in order to focus on the role of the network topology in power-grid synchronization: 
Hi = H, Di = D, P* i = P*, x' di = 0, \V*\ = 1 p.u. for all i, and all transmission lines have 
identical parameters. These assumptions together imply \E*\ = 1 p.u. for all i, = PP for 
generators, = —P^ for motors, KfJ^ = K for all pairs of nodes i and j that are connected 
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(and = 0 otherwise), and 7 ^ = 0. This corresponds to a version of the SM model for 

a homogeneous power grid in which all generators, motors, and transmission lines have identical 
parameters, but the network topology can be arbitrary. 


V. OBTAINING MODEL PARAMETERS REQUIRED FOR SIMULATIONS 

The parameters required to run a simulation of the synchronization dynamics using Eq. (2) are 
the system reference frequency Ur, the parameters of individual generators/motors, Hi, Di, and 
x' di , as well as the constants for the network dynamics equation, Ai, Kij , and 77 . The parameters 
ur, Hi, Di, and x' di must be given, while Ai, K^, and 7 ij are computed by solving the power flow 
equations (3) and (4) and performing Kron reduction. As described in Section II, a power flow 
solution is usually calculated given the parameters P* i and \V*\ at each non-reference generator, 
PJ i and Q\ - at each load node, and | V* | at the reference generator, as well as the admittance matrix 
Yo representing the physical network structure. The type of power system dataset most common in 
engineering is suitable for this calculation and thus contains the required parameters at each node, 
as well as sufficient information to compute Yo. Given such a dataset, standard power systems 
software, such as PST 20 and MATPOWER 21 , can be used to compute a power flow solution. The 
solution can then be used to build the matrix Yq incorporating the generators’ x' di , which can in 
turn be used to compute the constants A,f P , Kf?, and 7 fj* for the SP model. From Yq, P/ - , and 
Q\i, we build Y EN and then compute A EN , A EN , and 7 ™ for the EN model. For the SM model, we 
combine the motors’ x' d ■ with the matrix Yo to compute Y SM , which allows for the computation of 
A® M , KfJ^, and 7 ® M . The process of computing the model parameters required for simulations, as 
well as the dependency relations between all parameters involved, is illustrated in Fig. 2. We have 
implemented this process using MATPOWER, and our open-source software is freely downloadable 
from https://sourceforge.net/projects/pg-sync-models. The parameter computation and 
simulation of the dynamics can be performed by PST 20 for the EN model, but not for the SP or 
SM models. 

The dynamical parameters for generators and motors, Hi, Di, and x' di , are often not available. 
For unavailable values of x' di and Hi, we suggest estimating them using the strong correlation 
between each of these parameters and the active power output P* - of generator %, which was 
observed in Ref. 8 for the first three systems listed in Table I. (The values for these three systems 


20 



FIG. 2. Computation of the model parameters {Aj. K, l j , and 7 ij) given the input data for 
simulating network synchronization dynamics and the dependency relations among the input data 
and the computed parameters. The parameters required for the EN, SP, and SM models are 
enclosed by orange, blue, and green boxes, respectively. 

are available from Refs. 11, 24, and 38.) The parameters are estimated as x' di = 92.8(P* i )~ 13 and 
Hi = 0.04 Pgi, while imposing a maximum value of one for x' di and a minimum value of 0.1 for 
Hi, where P* t is measured in megawatts, x' di is in p.u., and Hi is in seconds. Our method of 
estimating x' d • is a refinement of a standard approach based on the observation that real values of 
x' d i tend to fall within a narrow range when expressed in p.u. with respect to the rated power of the 
generator/motor 11 ’ 19 . This observation, combined with the assumption that the rated power (which 
is unavailable in many datasets) is proportional to Ph, implies that x' di is inversely proportional 
to P*i when x' di is expressed in p.u. with respect to a reference power common to the entire 
network. In our previous work 8 we found that the available data for the first three systems in 
Table 1 follows our formula more closely than the standard approach. For Di, one may use 
Di = (D mt i + D et i + 1/ Ri)uji/P r = 50 p.u., a value that results from assuming that mechanical 
friction and the electrical effect of damper windings are negligible (D m j = D e) i = 0 ) compared to 
the effect of regulation by governors with Ri = 0.02 c cr/Pr = 0.02 p.u. 


21 






























TABLE I. Structural and dynamical properties of the systems considered. 


System 

Nodes 

Links 

Generators 

Loads 

Number of oscillators 

EN SP SM 

3-generator test system (Ref. 11) 

9 

9 

3 

6 

3 

12 

9 

10-generator test system (Ref. 24) 

39 

46 

10 

29 

10 

49 

39 

50-generator test system (Ref. 38) 

145 

422 

50 

95 

50 

195 

145 

Guatemala power grid a 

370 

392 

94 

276 

94 

464 

370 

Northern Italy power grid a 

678 

822 

170 

508 

170 

695 

678 

Poland power grid (Ref. 21) 

2383 

2886 

327 

2056 

327 

2710 

2383 


a Data provided by F. Milano (University of Castilla - La Mancha). 


VI. HETEROGENEITY OF MODEL PARAMETERS 


Given the input data for the system, the constants A,, Kjj. and 7 ^ for the network dynamics 
governed by Eq. (2) can generally be distributed heterogeneously across the network. Since Aj 
are related to the inherent frequencies to* through the same formula to* = to r(1 + Ai/ Di) for all 
three models, the Aj-heterogeneity quantifies the heterogeneity of the intrinsic properties of the 
generator dynamics. While the heterogeneity in the generator parameters is naturally a source of 
heterogeneity in A,, it is actually possible to have heterogeneous A, even in the complete absence 
of heterogeneity in the generator parameters Hi, Di, x' di , PL, and \V*\. The first two parameters 
actually do not enter into the calculation of A* (as described in the previous section) and thus do 
not affect the Aj-heterogeneity. The other three affect the values of Ai, but do not fully determine 
Ai, since the calculation also involves the admittance matrix Yo as well as P/, and Q/ i for the 
loads. In the following subsections, we first show that Aj-heterogeneity can indeed arise without 
heterogeneity of generator parameters using a small network example, and generalize it to larger 
networks. We then discuss the Aj-heterogeneity as well as the coupling matrices K tJ and 7 ^ for 
examples of realistic systems. 
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A. System with two generators and one load 


We consider the simplest case with which the heterogeneity of generator parameters can be 
discussed: a system with two generators, each connected to a single load by a transmission line, as 
shown in Fig. 1. For this system we have n g = 2, = 1, n = n g + ni = 3, and N = 2 n g + ni = 5. 

Assuming lossless and inductive transmission lines (i.e., r\ 3 = r23 = 0, X13 > 0, and X23 > 0), the 
power flow equations for this network read 
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( 22 ) 

(23) 

(24) 

(25) 

(26) 
(27) 


where B 0 11 = Im(Y 0 n) = V ~ Tfi > B022 = ^ and -^033 = -B011 + B 0 22- Setting 

0i = 0 and using it as the reference angle, we solve this set of five equations for five unknowns, 02, 
03, Q 1, Q 2 , | V31 to obtain a power flow solution. (The third equation does not involve unknown 
variables and simply imposes a constraint on the parameters.) The equations of motion for the 
three network dynamics models are: 


EN model: 
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SP model: 
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SM model: 
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Note that the phase shifts / yf p and 7^ M do not appear in the equations, since the assumption of 
lossless and inductive transmission lines implies that they are both zero, as discussed in Sec. IV B 
and Sec. IV C. 

Assuming that P*i, P* 2 , | V* |, and |V 2 *| are all nonzero, Eqs. (22) and (23) imply 

Agq 1^2^131 = sin(0i - fa) 

A* 2 |V?* 23 | sin (fc-fo) 

(regardless of the choice of the network dynamics model). Thus, even if the two generators’ output 
power and voltage magnitudes are identical, the steady-state phase angle differences (with respect 
to the phase of the load node) can be different if the parameters of the transmission lines are 
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different. This has an interesting implication for the EN model: the oscillators can have different 
parameters (i.e., A™ 7^ A 2 N , an d therefore different inherent frequencies) even when the generator 
they represent have identical parameters (i.e., P* x = P* 2 , Vj* = V 2 . and x' dl = x ' d2 ) if the 
transmission lines have different parameters. This is because we have A PN = P* — |E*| 2 G PN , the 
constant \E*\ depends on the steady-state phase angles, and G PN depends on the transmission line 
parameters. The effect is illustrated in Fig. 1C using the following concrete values of the system 
parameters: X13 = 0.04 p.u., X 23 = 0.047 p.u., 613 = 0.082 p.u., 623 = 0.098 p.u., P* 1 = P* 2 = 
1 p.u. = 100 MW, P£ :i = 2 p.u. = 200 MW, = 1 p.u. = 100 MVAR, V* = V 2 = 1 p.u. = 100 kV, 
x' dl = x' d2 = 0.1 p.u. (with all p.u. quantities with respect to the 100 MVA system base). If we 
lift the assumption of lossless transmission lines (purely imaginary impedances), the effect we 
have demonstrated for the EN model can also be observed for the SM model. Indeed, if we set 
ri3 = 0.002 p.u. and r23 = 0.007 p.u., while adjusting P* := P* = P 2 so that the power flow 
equation is satisfied and keeping all the other parameters the same, we obtain = 0.9217 and 
A 2 M = 0.8109 (with P* = 1.0055 p.u.) for the SM model. In contrast, the same effect cannot be 
observed for the SP model, since we have A? p = PP for all generators regardless of the network 
structure. 

The effect of network asymmetry demonstrated here for this small example system when using 
the EN and SM models suggests generalization to larger systems, as well as to other forms of net¬ 
work asymmetry, such as heterogeneities in the distribution of power demand, generator locations, 
and the connectivity of generators/loads in terms of the number of transmission lines. This will 
be discussed in the next section. 


B. Larger power-grid networks 

Consider a system with identical generator parameters: PP = P*, \ V* | = | V* |, and x' d i = x' d for 
alii = 1,... , rig . Consider the EN model. The constants Af N = P* — \E*\ 2 Gf^ are identical if and 
only if |E*| 2 G ,pn are identical. From Eq. (10), we see that requiring identical \E*\ constrains all Q*. 
Also, requiring G PN to be identical is equivalent to a global constraint through Eq. (14). We thus see 
that A PN will not be identical for a power flow solution in a generic setting. Nongeneric situations 
can arise if there is a symmetry in the system, such as the system with the star configuration 
in which all generators are connected to a single load through identical transmission lines. Even 
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with this star configuration, any heterogeneity in the parameters of the transmission lines will lead 
to heterogeneity in cj)* and Q* at the generator terminal nodes, and thus in \E*\ and A? N . A 
similar argument can be made for the SM model. We therefore conclude that a generic system has 
heterogeneous Ai under both the EN model and the SM model. 

We now study the heterogeneity of the model parameters for the set of realistic systems listed 
in Table I, which are used in Ref. 8 to study the stability of synchronous states against small 
perturbations. The data required for power flow calculations for these systems were obtained as 
follows. For the 10-generator system, known as the New England test system, the parameters were 
taken from Ref. 24. For the 3- and 50-generator systems, the parameters were taken from Refs. 11 
and 38, respectively. The data for the Guatemala and Northern Italy systems were provided by F. 
Milano (University of Castilla - La Mancha), and the data for the Poland system were provided as 
part of the MATPOWER software package 21 . The required dynamic data that are not available 
from the respective sources cited above were computed as described in Section V. 

Figure 3 shows the computed inherent frequencies ui* = wr(1 + Ai/Di) under the three models 
for each of the systems in Table I. Note that those realistic systems have heterogeneous generator 
parameters, which is a source of heterogeneity in A, in addition to the effect described above 
that results from the heterogeneity and asymmetry of the physical network structure. We see 
that the inherent frequencies can be significantly different from the system reference frequency 
cjr, with mostly higher frequencies for the generators and lower frequencies for the load nodes 
(except for the EN model, in which the load nodes are eliminated and thus have no well-defined 
inherent frequencies). This is consistent with the general picture that the voltage phase angles of 
the generators tend to rotate at higher frequencies and pull the frequency of the rest of the grid 
upward, while the load phase angles tend to rotate at lower frequencies and pull the rest of the 
grid downward. The variation of the inherent frequencies among generators and load nodes can be 
extremely large. Indeed, we find a generator with w* > 7 x 10 2 Hz for the 50-generator system, 
and a load node with uj* < —1.1 x 10 3 Hz (corresponding to a rotation in the opposite direction 
with frequency > 1.1 x 10 3 Hz) for the Northern Italy system. As noted in Sec. I, frequencies far 
from wr are not likely to be observed in reality, but these values of inherent frequency do reflect 
the inherent dynamical property of the oscillators representing generators or loads. The observed 
large variation in inherent frequency across the network may indicate a high level of frustration 
in the system in the sense that oscillators with vastly different frequencies are forced to rotate at 
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FIG. 3. The inherent frequencies oj* of the nodes in the systems listed in Table I under the 
EN (orange), SP (blue), and SM (green) models. The nodes are indexed so that i = 1 ,...,n g 
correspond to the generator internal nodes, the i = n g + 1,...,2 n g to the generator terminal 
nodes, and i = 2 n g + 1,...,JV(= 2 n g + ng) to the load nodes. The three types of nodes are 
separated by the two gray vertical lines at i = n g and i = 2 n g . The plots show u* of the generator 
internal nodes for all three models, the generator terminal nodes only for the SP model, and the 
load nodes only for the SP and SM models. The black horizontal lines indicate each system’s 
reference frequency. We note that the inherent frequency of some generators and loads are outside 
the range of the plots. 


a common frequency by the coupling. Comparing across different network dynamics models, we 
see similarity in the pattern of generator-to-generator variation, as well as noticeable differences in 
the inherent frequency of the same generator or load under different models. The difference across 
models does not always depend monotonically on the size of the system, as evident from comparing 
the generator frequencies in the 10-generator system, 50-generator system, and Guatemala power 
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FIG. 4. (A)-(C) Coupling matrix K = ( Kij) for the three different models of the 50-generator 

system (with n g = 50, ri£ = 95, n = 145, N = 195). (D)-(F) The corresponding phase shift matrix 
r = ( 7 ij). In each row, the common color scale used for all three matrices corresponding to the 
three models is shown on the right. The matrix elements are colored white if they do not affect 
the network dynamics (those elements for which i = j or K, : j = 0). 

grid. 

Figure 4 shows a visualization of the coupling matrix K = {Kij) and the phase shift matrix 
r = ( 7 for the 50-generator system. Recall first that the size of the coupling and phase shift 
matrices is different for different models: n g x n g for the EN model, N x N for the SP model (where 
N = 2 rig+ne), and nxn for the SM model (where n = n g +ne). Despite the fundamental differences 
in the dimensionality and in the type of equations used in the three models, we can identify some 
similarity within the coupling structure. We see in Fig. 4A that there is a non-trivial structure 
in the dynamical coupling among the generator internal nodes under the EN model. Although 
this coupling matrix is fully-populated, the strength of coupling varies widely and spans across 
many orders of magnitude, from 10 ~ 9 to order one. Notice that a similar coupling structure can be 














































identified in the middle box in Fig. 4B, corresponding to the sparse physical network connecting the 
generator terminal nodes in the SP model. We can also identify a similar structure in the top left 
box in Fig. 4C for the SM model. Most of the 7 ^ values are close to zero (red in Fig. 4D-F), which 
indicates that the coupling between oscillators i and j tends to make the oscillators synchronized 
with a small phase angle difference. A zero phase shift corresponds to an admittance (physical 
or effective, depending on the model) that has positive imaginary part and zero real part, which 
means that it represents a pure capacitive reactance. There are, however, some pairs of oscillators 
for which 7 ^ is significantly different from zero (green to blue in Fig. 4D-F). Such large phase 
shifts correspond to resistive or inductive components. 


VII. CONCLUSIONS 

We have presented first-principle derivations for three different network-level models of power- 
grid synchronization dynamics. These models are all based on the same fundamental equation of 
motion for mechanical rotation (the swing equation) and the same electric circuit representation 
(the classical model) of power generators and synchronous motors. The models are, however, 
clearly distinguished by the different sets of assumptions adopted for the loads. All three models 
are formulated in their full generality, which helps clarify how different choices of assumptions lead 
to different special cases of these models that have been used in previous studies. We have shown 
that even when we assume all the generators to be completely identical, the parameters of the 
phase oscillators representing them in the EN and SM models may be quite different. We have 
also shown that the oscillator parameters are vastly heterogeneous in a selection of test systems 
and real power-grid datasets. 

While the validity and appropriateness of the models depend on the purpose, the choice of addi¬ 
tional assumptions, and the system under consideration, these three models lay a solid foundation 
for the study of synchronization dynamics in complex power-grid networks. If the purpose is to 
isolate the effect of the network topology on synchronization, making simplifying assumptions (as 
has been done in many previous studies) to homogenize all other aspects of the system may be 
appropriate. If the purpose is to investigate the interplay between heterogeneities in the param¬ 
eters of the generators, loads, and transmission lines, it may be helpful to use a simple, random, 
and/or regular network topology. If the purpose is to understand the network structure and the 
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component parameters of realistic power systems, then keeping the details in a given dataset and 
comparing the synchronization properties with null models will provide the most useful insights. 
The models we considered here can serve as bases for stability analysis against small or large per¬ 
turbations, for which a variety of analytical and numerical techniques are available. As shown in 
Ref. 8 , stability against small perturbations can be fully addressed using the powerful conceptual 
framework of master stability functions 29 based on linear stability analysis and network spectral 
theory. In contrast, results on stability against large perturbations tend to be more limited in scope. 
If the relevant admittance matrix has zero real components, the problem can be tackled by energy 
function-based approaches 24 , which can also be used in combination with bi-stable representation 
of transmission lines for the modeling of the dynamics of cascading failure propagation in power 
grids 39 . Other approaches include estimating the region of attraction of a synchronous state 40 , 
quantifying the size of the region using a recently proposed measure called the basin stability 41 , 
and analyzing nonlinear modes using Koopman operators 42,43 . Such analyses of synchronization 
stability against small and large perturbations may help develop strategies for improving the exist¬ 
ing power grids by elucidating the relation between network structure and stability, which we have 
recently demonstrated to be highly non-intuitive in general models of coupled oscillator networks 44 . 

Understanding the dynamics of the models we considered here can potentially serve as a stepping 
stone for studying other types of instabilities relevant for power grids. For example, studying the 
so-called dynamic bifurcation problem 15 for these models, in which slow variation of the parameters 
A t . Kij, and/or 7 ^ drives the system toward a bifurcation point, can provide insights into the issue 
of voltage instability in power grids 46 . To capture the fast dynamics that follows the bifurcation, 
one may supplement the models we studied with additional equations that describe time-varying 
voltage magnitudes 47-49 . Such model-based studies of voltage instability complement model-free 
approaches for predicting the bifurcation point based on time series measurement 50 . Other causes 
for instability include fluctuations in power demand and failure of system components, and model¬ 
ing them as stochastic processes can help optimize response strategies for power grid operators 51 . 
With the anticipated transition to smart grids, there will be additional fluctuations in the system 
from new components such as intermittent energy sources and increased use of plug-in electric 
vehicles 52 . Such fluctuations will have impact on the synchronization stability of power grids, and 
its implications will be difficult to understand unless we understand the short-term dynamics of 
the models considered here. By providing comprehensive comparative analysis of these models at 
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the intersection of the physics of complex systems, network science, and power system engineering, 
we hope to contribute to a better understanding and control of the dynamics of existing power 
grids as well as to a better design of future power grids. 
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Appendix: Single generator connected to a large power grid 

Consider a generator connected through a single transmission line to a large power grid. Since 
the power grid is large, the dynamics of the generator does not affect the behavior of the rest of the 
grid much, and we thus approximate the voltage magnitude and the phase angular frequency at 
the other end of the line to be constant. This approximation corresponds to having the generator 
connected to a sort of “heat bath,” which can be modeled as an artificial generator that has 
infinite inertia and can absorb an arbitrary amount of active and reactive power without changing 
its terminal voltage magnitude and phase frequency (and thus is often called an infinite bus, where 
the term “bus” is used by power system engineers to refer to a network node). We measure 
the generator’s phase angle <5 relative to the phase of the voltage at the power-receiving end of 
the transmission line, which has magnitude |V)>o| and phase rotating at the reference frequency 
ujr. Combining the transient reactance and the transmission line impedance Z = R + jX to 
obtain a total admittance of Y = [R + j(X + a^)] -1 = |T^|e JQ! = G + jB (which defines a), 
the active power output P e from the generator at its internal node can be calculated as P e = 
\E*\ 2 G + \E*VooY\ sin(<5 — 7 ), where 7 := a — n/2. Equation ( 8 ) then becomes 

9 tj d 

- 6 H- 5 = A — A'sin(<5 — 7 ), (A.l) 

UR UR 

where A = P m — \E*\ 2 G and K = \E*V 00 Y\. Upon changing the variable to 6 := 5 — 7 and defining 
M := 2 H/ojr and D := D/lor, the equation becomes M9 + D0 = A — K sin#, which is the equation 
of motion for a forced damped pendulum, where 6 is the angle the pendulum arm makes with the 
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vertical, M is the mass, D is the damping coefficient, and A is the forcing torque. We thus see 
that a pair of stable and unstable equilibria exist if and only if A < K, and they are characterized 
by 6* = 0 and A = K sin 6* and correspond to the states of the generator synchronized to that of 
the rest of the grid at the frequency ujr. 

If D 2 > 4 M\jK 2 — A 2 , the effect of damping is strong enough, and 5 will converge exponentially 
to zero (corresponding to the generator’s instantaneous frequency converging to ujr) after a small 
perturbation is applied when the system is at the stable equilibrium. Otherwise, the damping is 
too weak, and the frequency will oscillate around uir with decaying amplitude with a characteristic 
frequency determined by the parameters. When D = 0, the equilibrium is neutrally stable and the 
(angular) frequency of oscillation around the equilibrium is ( K 2 — A 2 ) 1 / 4 / \[M , which is referred 
to as the “natural frequency” of the generator in the power systems engineering literature 11 . The 
natural frequency is different from the inherent frequency discussed in the main text, which gives 
uj* = lor + (P rn — | E*\ 2 G)/D (corresponding to <5* = A/D) and corresponds to the equilibrium 
frequency of the generator’s rotor in the absence of the coupling term (i.e., K = 0). The K = 0 
decoupling situation can in fact be realized briefly if a fault occurs on the transmission line near 
the other end and creates a short circuit to the ground. 

We emphasize here that, among the parameters that determine the dynamics of Eq. (A.l), the 
constants \E*\ and P m depend on how much power is flowing over the transmission line, in contrast 
to all the others, which represent the property of the system itself and are independent of the power 
flow. In a practical setting, the power flow over transmission lines is determined by solving the 
power flow equations (3) and (4), given the active power output and the voltage magnitude at 
the generator terminal nodes as well as the active and reactive power consumption at the load 
nodes. Given the voltage magnitude |V*| and the active power P* injected by the generator at the 
terminal into the transmission line in the single generator example of Eq. (A.l), Eqs. (3) and (4) 
(with n = 2) give 

Pg = cos(0 + a), Q*g = ~ lV J~ l sin(0 + a), (A.2) 

where cj) is the voltage phase angle at the terminal relative to that at the other end of the trans¬ 
mission line and Q* is the reactive power injected into the line. This can be used to express Q* 
in terms of the given parameters as Q* = y // |E*V' 0O | 2 /|.Z’| 2 — ( P* ) 2 . This, together with P m = P* 
and Eq. (10), determines the dependence of the parameters of Eq. (A.l) on the steady-state power 
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flow of the system, specifically as a function of P*, |V*|, and 
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